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Abstract: Tissue phantoms play a central role in validating biomedical 
imaging techniques. Here we employ a series of methods that aim to fully 
determine the optical properties, i.e., the refractive index «, absorption 
coefficient ju a , transport mean free path i , and scattering coefficient ju s of a 
T1O2 in gelatin phantom intended for use in optoacoustic imaging. For the 
determination of the key parameters jj a and I , we employ a variant of time 
of flight measurements, where fiber optodes are immersed into the phantom 
to minimize the influence of boundaries. The robustness of the method was 
verified with Monte Carlo simulations, where the experimentally obtained 
values served as input parameters for the simulations. The excellent 
agreement between simulations and experiments confirmed the reliability of 
the results. The parameters determined at 780 nm are n = 1.359(±0.002) , 
M ' s = 1 / f = 0.22(±0.02) mm 1 , jU a = 0.0053(+0.0006-0.0003) mm 1 , and 
fi s = 2.86(±0.04) mm ". The asymmetry parameter g obtained from the 
parameters £ and ju' s is 0.93, which indicates that the scattering entities 
are not bare Ti0 2 particles but large sparse clusters. The interaction between 
the scattering particles and the gelatin matrix should be taken into account 
when developing such phantoms. 
© 2012 Optical Society of America 

OCIS codes: (170.3660) Light propagation in tissues; (170.6935) Tissue characterization; 
(290.1990) Diffusion; (290.4210) Multiple scattering. 
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1. Introduction 

The development of diagnostic imaging techniques commonly entails the use of "phantoms", 
which try to mimic human tissue as closely as possible. They are conceived in such a way that 
their optical properties are similar to those of tissues and, therefore, constitute a precious tool 
which provides a well-defined environment to study light propagation. Quoting a recent 
comprehensive review by Pogue and Patterson [1], these phantoms "are used for a number of 
purposes, including: (i) initially testing system designs, (ii) optimizing signal to noise in 
existing systems, (iii) performing routine quality control, (iv) comparing performance 
between systems." 

Our group is conducting research in optoacoustic (OA) imaging [2] where phantoms are 
widely used based on hydrogels (such as gelatin, agar, polyvinylalcohol) or PVC plastisol and 
other resins (see [1,3-5] and references therein). We opted for gelatin phantoms because they 
are relatively easy to produce and have the advantage of being adjustable to various purposes: 
(i) the phantom's optical properties can be customized by incorporating scattering agents (like 
Ti0 2 or intralipid) or absorbing elements (like India ink, or any dye of interest), and (ii) the 
phantom's mechanical/elastic properties can be modified by varying the fraction of gelatin. 
Elastic properties especially play a valuable role: indeed, the combination of OA imaging with 
ultrasound elastography holds promise for reduced clutter and increased imaging depth in the 
framework of clinical imaging [6,7]. 

Our gelatin phantoms are produced in a two-step approach. First, a phantom precursor 
stock is prepared by dissolving gelatin in water and by incorporating titanium dioxide (Ti0 2 ) 
as the scattering component. During the second step, different phantoms can be developed 
from this precursor stock by adding various amounts of supplementary scattering agents, 
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absorbing structures, and acoustic scatterers (such as flour) to obtain more realistic optical and 
acoustic properties for combined OA and ultrasound imaging. 

Accurate knowledge of the optical properties of the phantom precursor stock is essential 
within this process. Determining the optical properties of such phantoms or tissues is not an 
obvious task and has been subject to intensive research over the past decades (see e.g. the 
review by Kim and Wilson [8]). In this paper, we employ a series of methods by which we 
were able to fully characterize the optical properties of the precursor stock (which is simply 
referred to as "phantom" or "sample" throughout this manuscript) at a wavelength of 780 nm. 
This wavelength is particularly appropriate since it lies within the diagnostic window where, 
owing to low optical absorption, OA imaging can achieve several centimeters of imaging 
depth [9]. 

Experiments 
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Fig. 1. Diagram summarizing the different stages of the work presented in this paper: the 
optical properties (refractive index, transport mean free path, absorption coefficient, scattering 
coefficient, and scattering law) of the phantom were experimentally determined by the 
indicated set of experiments. Monte Carlo simulations were performed ex-post to validate the 
experimental results. 

The way this paper is organized is outlined in the diagram of Fig. 1. After describing how 
the gelatin phantom was produced, we describe how the optical properties of the phantom 
were determined at 780 nm, as follows: 

a) its refractive index n was determined with an Abbe refractometer; 

b) its transport mean free path i* and absorption length l u =l///„ (where /j a is the 

absorption coefficient) were measured by performing time of flight (TOF) 
measurements within the phantom; 

c) its scattering coefficient /j s was assessed from extinction measurements with an 

absorption spectrometer; 

d) and its scattering law (which models the way the phantom's scatterers spatially 

redistribute the light) together with its corresponding asymmetry (or anisotropy) 
parameter g were evaluated from goniometric measurements. 

In the last part we present Monte Carlo (MC) simulations performed to validate the 
experimental results: we simulated the TOF experiments by providing the experimentally 
determined optical properties of the phantom as inputs to the MC program. 

2. Gelatin phantom preparation 

2.1. Ingredients and preparation of the phantom 

The phantom precursor mixture consists of 16.7% weight fraction of gelatin dissolved in 
water and 1 gram per liter (g L -1 ) of Ti0 2 for optical scattering. The Ti0 2 concentration was 
chosen according to the data in [4], which suggest a reduced scattering coefficient of 0.4 mm -1 
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at 1 g L -1 . This is below the range 0.5-1.3 mm -1 typically found in human breast tissue 
[10,11], but we anticipate the need for further additives, including e.g., flour for acoustic 
contrast, which will modify the scattering properties. 

The Ti0 2 particles were obtained from Millenium Chemicals (type Tiona RL11A). They 
consist of 99% rutile, with a refractive index of 2.4884 at 780 nm [12]. Their size distribution 
lies between 180 nm and 400 nm diameter, with a clear peak at 250 nm (weight percentage, 
determined by disc centrifugation). 

The phantom in its precursor form is prepared as follows: 

• With the help of an agitator, 1 kg of gelatin is mixed with 5 L of a stock solution 

prepared with deionized water and composed of 1 g L -1 Ti0 2 and 1 g L -1 
methylparabene (which serves as preservative), so that the final phantom contains 
16.7% weight percentage of gelatin. (For the dilution series used in Sections 4 and 5, 
the aqueous Ti0 2 dispersions are pre-diluted to the required Ti0 2 concentrations.) 

• The mixture is heated to 60°C, so that the gelatin gets dissolved and the mixture 

becomes homogeneous. 

• After heating the mixture up to 60°C, it is left to cool down to 50°C, with the agitator 

running throughout the whole heating and cooling processes, in order to maintain the 
homogeneity of the mixture. 

• Finally, the resulting bulk solution is poured into a glass aquarium (of size 

1 80 mm x 240 mm x 2 10 mm : see Fig. 2) where it cools down to room temperature 

and becomes solid. As explained in Section 3, the phantom is prepared in such a 
large container in order to avoid boundary condition problems during the TOF 
measurements 

2.2. Refractive index 

For the evaluation of TOF measurements (Section 3) and for the MC simulations (Section 6), 
it is necessary to know the refractive index n of the phantom at the wavelength 
of X = 780 nm . After having prepared a phantom sample, we measured at the room 
temperature of 22°C both its refractive index n D « 1 .363 at the sodium d-line wavelength 
X D = 589.29 nm and its corresponding Abbe number v = II (n F -n c ) = 56.72, where n F and 
n c are the refractive indices at the wavelengths X F = 486.13 nm and X c = 656.27 nm . These 
two values were then used to determine the two parameters A and B of the two-term form of 
the Cauchy extrapolation formula n{X) = A + B I X 2 . This law gives the refractive index of the 
phantom at 780 nm: n = 1.359 ±0.002 . The measurements were done with an Abbe 
refractometer (Type IT by ATAGO). Since the refractrometer is calibrated with water, we use 
the error estimate of the IAPWS formulation for pure water in the wavelength range 700-1 100 
nm[13]. 

3. Determination of the transport mean free path and absorption length 

3.1. Principle of the measurements 

Our technique is a variant of position and time resolved diffuse reflectometry [14], but we 
place photon injection and detection deep inside the scattering medium, in order to minimize 
the influence of the boundaries [15-17]. The principle is illustrated in Fig. 2. Photons 
generated by a pulsed picosecond diode laser (LDH-P-785, PicoQuant) are delivered into the 
phantom through a multimode optical fiber (50 mm core diameter) and injected at a point S 
well inside the scattering medium. A second fiber collects the scattered light at a point D at a 
distance r from S and guides it to a single photon detector (idl00-50, ID Quantique). The time 
of flight (TOF) of the individual photons from S to D is measured by a fast timer (SPC 530 
PicoQuant) triggered by the laser pulse. The resulting TOF distribution is characteristic of the 
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optical properties of the material through which the light propagates. Thus, by fitting an 
appropriate analytical model to the recorded distribution curve, the optical properties of the 
phantom can be extracted. 
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Fig. 2. Schematic representation of the principle of the time of flight (TOF) experiments: the 
light source - a picosecond laser - injects the light into the phantom at a point S; the individual 
photons are then subject to multiple scattering in the phantom (random photon paths are 
represented in red) until they exit the phantom or get absorbed in it. The photons that reach the 
detection point D (placed at a distance r from S) are detected and their individual TOF from S 
to D recorded in a histogram. This experiment is reproduced for different source-detector 
separations r and the obtained TOF distributions are then fitted to the diffusion model. 



3.2. Diffusion model 

We employ the photon diffusion model [18], which is an approximation to describe the light 
transport in strongly scattering media (i.e., where ju t » jU a ) and over large distances such 
that r » f . According to this approximation, the probability P(r, t)dr 3 dt to detect a photon 
in a volume element dr 3 and time interval dt obeys the diffusion equation: 



dP(r,t) 

dt 



■ DV 2 P(r, t) = -aP(r, t) + S(r, t) 



(1) 



where D p is the photon diffusion coefficient and a p is the absorption rate. D p and a p are 
related to the phantom's transport mean free path f and absorption length t respectively: 



and a p =- 



(2) 



Here c, is the light speed in the medium. (In the biomedical photonics community one often 
uses the reduced scattering coefficient fj! s =\lfw&fi a =\lt a . Note that we assume 
jU s » jU a .) The Green's solution of Eq. (1) is 
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This simple solution applies in infinite space and for a 5-source in space and time: 
S(r,t) = S(r)S(t) . To approximate the former condition, we employ a large sample container, 
so that photons are injected and detected far from boundaries. To regard 5(r) as a point 
source, we require that ds <s t , where ds is the characteristic size of the source. This 
condition is easily met by using an optical fiber as an illuminator. 

The time dependence of the measured distribution M(t) results from the convolution of 
P(f) with the instrument response function IRF(f) that includes the finite width and shape of 
the laser pulse and the timing jitter of the detection: M = IRF ® P . To determine the least 
square estimates of the parameters D p and a p , we apply the standard iterative reconvolution 
technique, using the Levenberg - Marquardt algorithm provided by Matlab. A third parameter 
is needed to account for the transition from the ballistic photon transport to diffusion that 
occurs within approximately a distance of 21* from the injection point. We include in the 
fitting procedure a (manually) variable shift t d of the measured TOF-distribution with respect 
to the measured IRF. This shift does not play a significant role if r » f . Since f is not 
known a priori, we must carry out a series of measurements at varying source-detector 
separations. These estimated values should converge to the "true" values as r increases, so 
that most of the photon paths are long enough to assure the validity of the diffusion 
approximation. 

3.3. Setup and measurement procedure 

The phantom is contained in a soda-lime glass aquarium ( 1 80 mm x 240 mm x 2 1 0 mm ). The 

setup was designed to perform measurements at different source-detector separations without 
deforming either the phantom or the fibers (see Fig. 3). The multimode fibers (core diameter 
of 50 mm), both on the laser and the detector side, are guided through hypodermic needles 
(diameter 0.8 mm) in order to reach the inner parts of the gelatin phantom with minimal 
damage and better accuracy. Up to eight such optodes can be positioned in the orifices of a 
holder-plate and fixed perpendicularly to the phantom's top surface. The heads of the needles 
are placed on the holder so that their tips are deep inside the phantom. This allows both the 
injection and the collection of the light far from the phantom-air boundary (typically at a 
depth of 40 mm beneath the top surface). To avoid distortions caused by the insertion of the 
needles, the phantom was poured in its liquid state into the aquarium where the holder and 
fibers had been positioned beforehand; the phantom then solidified with the needles already in 
place. The accuracy of the positioning of the fiber tips was inspected by taking photographs 
prior and after the solidification of the gelatin: we found a maximum deviation of 0.6 mm 
from the intended separation. 

Each pair of fiber/orifice on the holder plate corresponds to a different source-detector 
separation r. Thus, measurements at different source-detector separations are carried out by 
successively selecting different pairs of fibers and connecting each fiber to the laser and the 
single photon counting module. Prior to each measurement the IRF corresponding to the 
chosen pair of optodes was determined by measuring the time-resolved scattering from a sheet 
of white paper immersed in water. With the present laser and photon detector the width of the 
IRF is small when compared with the width of the TOF distribution Mit). Therefore, the 
precise shape of the IRF does not play a crucial role. What matters, however, is the timing, 
i.e., the shift between IRF(t) and M(t). Strictly speaking, for a precise determination of the 
IRF, one should place the fiber tips close together so that the photon path length is shorter 
than if . However, this exact positioning is not easy to achieve (for example, we must insert 
filters to adjust the photon-count rate), and, therefore, in the data evaluation we shall rely on 
the adjustable parameter r d to compensate the resulting inaccuracy of timing. An example of 
the best fit of the measured data by the diffusion model is shown in Fig. 4. 
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Fig. 3. Detail of the setup used for the TOF measurements, (a) The gelatin phantom (5 L) is 
contained in a soda-lime glass aquarium, (b) and (c) In the holder-plate, the optodes can be 
placed at eight positions. Each pair of optode/orifice on the holder-plate corresponds to a 
different source-detector separation r. The tips of the needles are inside the phantom, at a depth 
of 40 mm beneath the phantom/air interface. 



Diffusion model 

MC simulation convoluted with experimental IRF 




Time [ns] 

Fig. 4. Example of the fit of the measured data (thin solid line in black) with the diffusion 
model (thicker line in red). The corresponding IRF is shown in green dots. For comparison, we 
also include the fit with the corresponding MC simulations (see Section 6) with the same 
parameters i and I (thicker line in blue). 

3.4. Results and discussion 

In Fig. 5 we compile the results of five series of measurements on five different phantoms. 
First of all we observe that £ is rather large, larger than the transport mean free path in 
human breast (typically around 1 mm [10,11]). Correspondingly large is the influence of 
inaccuracies and limitations of the diffusion approximation, which results in the data scatter at 
small distances r. From an analysis of the influence of primary experimental errors on the 
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least square estimates of the parameters I and t a we found that the overall variance is 
composed of three contributions: 

i) Accuracy of the fiber positioning, whose influence on t depends on the optode 

distance r. For r = 20 mm , an error in r of 1 mm results in 0.4 mm error in the 
estimation of I , at r = 50 mm , the resulting error is only about 0.15 mm. The 
impact of the positioning error on ^ a is negligible. 

ii) Accuracy in the determination of the parameter r d i.e., the shift between the IRF(t) 

and M(t). With an optode distance r = 20 mm , an error in v d of 50 ps will lead to 
0.5 mm misestimation of i , at r = 50mm the resulting error in i decreases to 
about 0.3 mm. In addition, the parameter r d accounts for about 50% of the variation 
of l a . 

iii) The remaining variation of l a and C (see e.g. the outlier, phantom 3 at 43 mm 
distance) is probably due to imperfect reproducibility of the fabrication and, in 
particular, to the inhomogeneity of the phantom. One likely cause of the 
inhomogeneity is the sedimentation of the Ti0 2 particles during the 
cooling/solidification process. 
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Fig. 5. The values of the transport mean free path (top) and absorption length (bottom) obtained 
from the TOF measurements with five different phantoms. The vertical dashed line is a "guide 
to the eye" indicating roughly the validity limit of the diffusion approximation. 
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Nevertheless, the estimates of I of the five phantoms converge with a reasonable consistency 
to the values t =4.5 ±0.5 mm and t a =190 + 15 mm, as r increases beyond approximately 
six times f . 

The transport mean free path of almost 5 mm (reduced scattering coefficient ju' s of only 
0.2 mm -1 ) is substantially larger than the typical values reported for elastomer based 
phantoms with the same Ti0 2 concentrations (DMSO 0.9 mm [19], PVCP 2.5 mm [4], 
Polyurethane 1.25 mm [20]). On the other hand, much larger values for t have been reported 
for a comparable gelatin based system: Pogue et al. [21] measured ju[(c) of Ti02 in 7% 
gelatin (at A = 690 nm ) and found a linear concentration dependence ju' s (c) as 0.14 -c ( ju[ in 
mm -1 , c in g L -1 ) up to c = 40 g L -1 . Their data suggest that f is as large as 7 mm for 1 g L -1 
Ti0 2 (unfortunately, the particle size is not specified). 

4. Scattering coefficient and asymmetry 

4.1. Extinction measurements 

The extinction coefficient /j e of the phantom was determined at 780 nm using a Perkin Elmer 
LAMBDA™ 750 absorption spectrometer. With the absorption coefficient ju a =l/( a as 
small as 0.005 mm -1 , we expect that most of the extinction is due to scattering, i.e., ju e = /j s . 
If the concentration of the TiC>2 solution is low enough that multiple scattering can be 
neglected, then the transmission follows Beer-Lambert's law / = I 0 exp[-d jU s (c)] , where d is 
the thickness of the sample, 7 0 is the intensity of the illuminating light and / the intensity of 
the transmitted light. Furthermore, ju s (c) depends on the concentration c of the Ti0 2 particles. 
In the simplest case, i.e., when concentration dependent clustering of titania particles can be 
neglected, ju s is proportional to the TiC>2 concentration c in the sample: 

ju s (c)=ca s Ti02 (4) 

where <j s TjQ1 is the scattering cross-section. Thus, the absorbance A = -log 10 (7 / /„) is 
expected to be linear in d and c: A = log I0 (e) <J s cd . The experimental data shown in Fig. 6 
show that A is indeed linear in both, the sample thickness d and concentration c, which 
indicates that we indeed measured in the single scattering regime and that the cross-section of 
the individual scattering entities does not change with the concentration. Moreover, the linear 
regime characteristic of the Beer-Lambert's law extends up to the concentration c = 1 g L~ , 
which we used to prepare the phantom. 

The scattering coefficient at a TiC>2 concentration of 1 g L -1 as obtained from the linear fit 
to the data in Fig. 6 is ju s (c = 1 g L~') = 2.86 ± 0.04 mm 1 . The error estimate represents the 
standard statistical error of the least square parameter estimation. 




TiO, concentration [g-L '] TiO, concentration [g-L '] 

Fig. 6. The absorbance A of TiO^ solutions (containing 16.7% weight percentage of gelatin) as 
the function of concentration c (left: sample thickness d = 0.1mm; right: d = 1mm). 
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4.2. Asymmetry parameter g 



Knowing the scattering coefficient fj. s , we can determine the phantom's asymmetry parameter 
g = ^cos(60) , i.e., the average of the cosine of the scattering angle 9 , which is related to the 
transport mean free path and the scattering coefficient according to 



1 



1 



(5) 



K (1 

Using f = 4.5 + 0.5 mm, and /u t = 2.86 + 0.04 mm" 1 , we obtain g = 0.93 ±0.01 . This 
large asymmetry indicates substantial forward scattering, which is much more pronounced 
than one would expect from Mie scattering on Ti0 2 particles. The largest g one can get with 
isotropic spherical particles of the given refractive index ratio would be about 0.7 (regardless 
of the size of the particles!). Moreover, finite element calculations of Thiele and French [22] 
indicate that the effect of non-spherical shape and anisotropy is not sufficiently large to 
explain the large asymmetry. Apparently, the scattering entities are not bare isolated Ti0 2 
particles. We speculate that the effective scatterers are some large (though sparse) Ti0 2 - 
gelatin clusters formed by micro-phase separation during the solidification of the system. 

5. Assessing the scattering law 

The ultimate characterization of a phantom tissue would consist of the complete determination 
of the scattering law p(6,(f>) (also called "phase function") characterizing the angular 
distribution of scattered photons in a single scattering event. With samples exhibiting strongly 
asymmetric, nearly forward, scattering, this is a rather difficult task, which requires 
specialized instrumentation. Here we report preliminary results obtained with a commercial 
goniometer originally intended for dynamic light scattering measurement (ALV GmbH 
Langen). The available hardware does not allow a precise determination of p{9) at small 
angles, but is sufficient to check the consistency of the results obtained so far. 

5.7. Goniometric measurements 

The geometry of the scattering experiment is sketched in Fig. 7. The light source was a 780 
nm CW diode laser (Laser Components). The scattered photons were collected with a 
multimode fiber receiver focused at the center of the sample and counted with a single photon 
detector. The sample is immersed in a toluene filled cylindrical vat, to minimize geometrical 



Laser 
780 nm 



Single- Mode 



780 nm >»• vFiber 

\y Polarizer 



( 1 ) : Surrounding cylinder filled with toluene 

o 8.5 cm 

(2) : Internal cylinder containing 

the TKVgelatin sample 
o 2 cm 




Photo-multiplier 



Fig. 7. Top view of the light scattering experiment. The light emitted by the laser is vertically 
polarized. We measured the light scattered by Ti0 2 suspensions (containing 16.7% weight 
percentage of gelatin) of different concentrations at various angles 40°<#<140°. 
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errors resulting from refractive index mismatch. A polarizer-analyzer pair selects the SS- 
component (perpendicular to the scattering plane) of the phase function. 

An accurate determination of the single scattering law of the strongly scattering sample 
would require an extrapolation to a sample cell whose dimensions are smaller than 1/ ju t , i.e., 
smaller than 0.4 mm. This is not practicable with the given setup and even less with a 
heterogeneous sample. However, the observed linearity of ju s (c) and ju[(c) (see [21]) 

indicates that the structure of the scattering entities does not depend on the Ti0 2 
concentration. This suggests a simpler alternative: We carried out measurements on Ti0 2 
dispersions (each containing 16.7% gelatin) with decreasing concentrations, and from this 
series we extrapolated the trend of the scattering law. Examples of the extrapolation procedure 
are indicated in the inset in Fig. 8. 

For each Ti0 2 -gelatin dispersion, the scattered intensity / meas (#) was measured at 

different angles 40° < 9 < 140° . (The apparatus is not optimized for forward scattering and 
therefore we restrict the measurements to scattering angles larger than 40°). The measured 
signal / meas (#) is normalized, so that I n (9) is obtained as follows: 

j /Q\ _ Aneas (^) ~ ^blank (^j (f\ 

A„,(#) 

Here I Uwk (9) is the blank signal, measured with water in place of the Ti0 2 -gelatin sample 
and I tol (9) is the reference signal measured with toluene. Since 1(9) oc ju s p(0) , the 
normalized signal can be expressed as: 

P„A#) 

ju s . lol and p tol (9) are the scattering coefficient and the scattering law of the toluene, 
respectively and jj s (c) is the concentration-dependent scattering coefficient of the phantom as 
previously determined. Toluene scatters as an almost perfectly isotropic dipole scatterer 
(Rayleigh scattering), thus: p tol (9) = 3/8;r . (The theoretical background of polarized Rayleigh 
scattering is discussed for example in [23].) Thus, the scattering law of the Ti0 2 suspension 
can be determined according to 

p(9\c) = ^-^-I n (9) (8) 

%K jU s (c) 

Moreover, jU t lol is identical to the so-called Rayleigh ratio R tol of the toluene. From the 
data provided by Moreels et al. [24], we determine the Rayleigh extrapolation formula 
tf, o ,(/L) = 1.65xl0 8 //l 4 (R tol in m" 1 and X in nm), which leads to n s tol =4.46xl0~ 4 m 1 . In 
Fig. 8 the extrapolated values are indicated with symbols "+". Note, that since jU s tol , p tol (9) 
and jU s are known, an adequate model should match not only the shape but also the 
magnitude of the experimental curve without adjusting a normalization pre-factor. 

5.2. Modeling of the scattering law 

As the first guess, an obvious model would be the Mie scattering law corresponding to nearly 
spherical Ti0 2 particles with a certain size distribution. However, as already mentioned, 
disregarding the size of the particles, the largest value of g that can be obtained with Ti0 2 
dispersed in the gelatin gel is 0.7. This is much smaller than the previously determined g = 
0.93. Thus, we need another scattering law to model the proposed clustering of the Ti0 2 in the 
gelatin. To account for the interference between the individual scatterers forming a cluster, we 
resort to the Rayleigh-Debye approximation. (Which is also suited to describe light scattering 
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in tissue in a more realistic fashion than Mie scattering. See Ricka et al. [23].) In the 
Rayleigh-Debye approximation one combines polarized scattering from elementary Rayleigh 
scatterers with a structure factor that characterizes their mutual arrangement. We employ the 
polarized version of the generalized Henyey-Greenstein scattering law (pvgHG) (discussed in 
more details in Refs. [23,25]), which assumes the following "fractal" structure factor F: 



F(q | L,v) °c- 



1 



where v > — 1 



(9) 



[l + CL|q|)T 

Here q is the scattering vector, |q| 2 = 2k 2 [l-cos(6?)] , where k is the wavenumber in the 
medium and 9 the scattering angle. The pvgHG structure factor has two modifiable 
parameters: the correlation length L characterizes the extent of correlation between the 
elementary scatterers, i.e., roughly the size of the proposed clusters. The power law exponent 
v is related to the shape of the correlation, i.e., roughly the shape of the clusters. The case 
v = 0 corresponds to the well-known Henyey-Greenstein phase function widely used in the 
field of tissue scattering. (Note, however, that we also include polarization. For more on these 
topics see Section 7.3 in [23].) The correlation length L is not yet known, but we already have 
an estimate of the asymmetry parameter g. Therefore, for the convenience of the data analysis, 
we numerically convert the parameters {L,v} into 



. (1) pvgHG {#=0.90, v=0) 
(2) pvgHG [g=0.93, v=0( 



. (3) pvgHG tg=0.93, i>=-0.3) 



0.05 



0.04 



■Z_ 0.03 

CD 



0.02 



0.01 




0 0.005 0.010 0.01 5 0.020 0.025 

TiO, concentration 



[gram per liter] 



0 SO 100 

Scattering angle 0 |degrees] 

Fig. 8. Estimating the scattering law of the Ti0 2 scatterers immersed in gelatin: the pvgHG 
scattering law model is fitted to the extrapolated data (represented by "+" symbols in the main 
figure) corresponding to infinitesimal concentrations of T1O2, i.e., where the single scattering 
regime applies. The fitting is achieved by modifying the set of parameters {g, v}. The 
extrapolation procedure which gives the values corresponding to infinitesimal concentrations is 
indicated in the inset figure. 

Thus, the pvgHG model is fitted to the experimental data by manually modifying the set of 
parameters {g,v} . (Recall that{g,v} are the only adjustable parameters; there is no adjustable 
pre-factor.) The examples in Fig. 8 indicate the sensitivity to the parameter g. The model with 
g = 0.90 and v = 0 (curve 1) fails to give an accurate description of the data but with g = 0.93 
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(curve 2) is more satisfactory. A further improvement is achieved by varying the parameter v. 
A reasonably good fit is obtained with [g = 0.93, v = -0.3} (curve 3). 

From g = 0.93 we estimate a correlation length L in the order of 10 |im. The exact value 
depends on the parameter v, but the data quality does not allow its unambiguous 
determination. Nevertheless, the scattering law is consistent with the previous results from 
Sections 3 and 4. 

6. MC simulations 

In the last stage of the investigation, MC simulations were performed to assess the robustness 
of the experimental results. We employ a simulation tool that we recently developed [25] to 
study polarized light propagation in complex geometries. While polarization is not important 
in the diffusion limit, we shall take advantage of the following features of our software: 

i) A user interface allows creation of complex shapes in the simulation space through 
which the ray tracing algorithm follows the paths of photons across the interfaces and through 
the different regions, ii) The user can choose the optical properties for different regions of the 
sample, including the interface properties (e.g. Fresnel reflections) and scattering laws. 
(Currently Mie and pvgHG are implemented.) These features allow transposing the TOF 
experiment into a realistic simulation (compare Fig. 9 and Fig. 2). The optical properties of 
the phantom and its container that were given as inputs to the simulations are listed in Table 1 . 

To simplify the simulation design and to speed up the computation we permit ourselves 
certain simplifications of photon launching and detection, which should not affect the results 
in the diffusion limit. Specifically, we neglect the presence of the needles, because they are 
made of non-absorbing material and thus create only a minor disturbance. 

Photons are simply launched downwards at a position corresponding to the tip of the 
source fiber, 40mm beneath the phantom's top surface. The individual photons are launched 
with randomly chosen initial polarization (the so-called unpolarized light) to mimic the 
depolarization of the laser light occurring in a multimode fiber. 




Number of 
absorbed/detected photons 



Light injection/detection 

at 40mm beneath 
the phantom/air interface 



Fig. 9. Screenshot from the MC simulation: the simulation tool was used to reproduce the TOF 
measurements presented in Section 3 (depicted in Fig. 2). The colored disc visualizes the 
number of photons detected by the individual detector rings. 
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The detection is simulated as a concentric array of detectors (r resolution of 1mm), which 
register incoming photons without absorbing them (see Fig. 9). We used two versions of this 
detection scheme. In the first version, a photon is registered when it crosses the horizontal 
plane of the receiving fiber tip within a numerical aperture NA = 0.87. The fibers used during 
the measurements have NA = 0.2, but increasing the NA reduces the simulation time. The 
second detection scheme registers a photon when a scattering event occurs within a narrow 
slab (thickness 1 mm) adjacent to the detector plane, again into the aperture of 0.87. (To 
compare the resulting TOF distributions for different scattering laws, they must be normalized 
by dividing with fi s .) This second detection scheme provides approximately a six fold 
improvement of counting statistics but is only appropriate for distances within the validity of 
the diffusion approximation. 

All the simulations presented in this paper have been performed on a Linux PC equipped 
with an Intel Pentium D computer processor with a clock speed of 3.40 GHz and a cache size 
of 2048 KB. The path generation routine written in C + + has been compiled with a GNU 
compiler (gcc version 4.3.2). 



Table 1. The optical properties (at X = 780 nm) assigned to the three different materials 
used in the experiment: the quantities listed here correspond to the experimental results 
and were provided as input parameters to the MC simulation program 



Parameters 


Phantom 


Glass Aquarium 


Air 


Refractive index 


1.359 


1.52 


1.0 


Transport mean free path t inpul 


4.5 mm 


00 


00 


Absorption length £ a i „ put 


190 mm 


oo (no absorption) 


co (no absorption) 


Scattering coefficient 


2.86 mirf 1 


co (no scattering) 


co (no scattering) 



6.1. Testing the diffusion limit, influence of scattering law 

In a first series of simulations we tested whether the setup indeed fulfills the conditions for the 
applicability of the diffusion approximation. The results for two optode distances r are shown 
in Fig. 10. The transport mean free path and absorption length are set to the previously 
determined values f input = 4.5 mm , l a input = 190 mm , but we vary the scattering law (see inset 
in the bottom of Fig. 10) and its parameters. Obviously, at r = 5mm the diffusion 
approximation fails to describe the simulated data (which is not surprising, see e.g [26].), but 
for 36 mm distance (%(*) we are already well within the diffusion regime. Moreover, at 
r - 36 mm there is no way to distinguish the scattering law used in the simulation, since t 
and l a are the only parameters that matter in the diffusion limit. However, the influence of 
the scattering law becomes apparent at an optode separation of about 1 1 , despite the fact that 
all simulations are done with the same £ and l a . Apparently, the asymmetry g is important 
on the length scale f and on the time scale that corresponds to a mean photon path length of 
1 1 . The asymmetry g as well as the form of the scattering law affect angularly resolved 
reflectance experiments [27,28]. 

6.2. Analyzing the simulated data 

In Fig. 4 the outcome of the simulation (pvgHG law { g = 0.93, v = -0.3 }) is compared with 
the measured data obtained at the optode separation r = 36 mm . The simulation is convoluted 
with the IRF of the experiment, i.e., we regard the simulation as a model function M(t), an 
alternative to the diffusion model. Both models match the measurement within the statistical 
error. However, there is a slight mismatch between the diffusion model and the simulation 
model with the same parameters. This could indicate that at r = 36 mm we are not yet within 
the range of the validity of the diffusion approximation. 
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Time [ns] 

Mie [g-0.62, n v Jn nmt J pvgHG fe=0.62, v=Q) 

Diffusion model pvgHG lg=0.93, i'=-().3} 



10 




Time [ns] 

Fig. 10. Examples of the MC simulations where we varied the scattering law and its 
parameters. The curves in the inset of the bottom figure are polarization averaged. The data 
were generated in an overnight simulation including 60 million photons for each scattering law 
(total number of scattering events: 47 x 10 9 ). The simulated curves show the raw, non- 
normalized and unsmoothed data with a time-resolution of 2 ps. At small source-detector 
separations (top), the limitation of the diffusion model (red line) is clearly visible but the 
differences between the scattering laws are hardly discernable in the semi-log plot. These 
differences are significant at short times, as shown in the linear plot in the inset of the top 
figure. (The small bump at 25 ps is probably a commensurability artifact due to the finite 
resolution of the ring detector.) At large source-detector separations (bottom), the applicability 
of the diffusion model becomes apparent. At both distances, the diffusion model is shifted by r d 
=50ps. 

To further elucidate this issue, we analyzed the simulation data by fitting the diffusion 
model. Note that thereby we regard the simulation data as equivalent to a measurement, the 
only difference from a real measurement is that we have now a perfect estimate of the IRF, 
namely a 8-peak positioned at t = 0. The fit parameters are the transport mean free path f , 
absorption length i a and the shift x d that accounts for the transition from ballistic to 
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diffusive regime of photon propagation. Examples of the results for four optode distances are 
summarized in Table 2. The estimates of the parameter uncertainty are derived from visual 
inspection of the quality of the fit. (Statistical analysis in terms of % 2 would not be appropriate, 
since we expect a systematic discrepancy between the diffusion model and the data for short 
photon paths.) 

The optimal shift r d = 50 ps , corresponding to 1 1 mm of the photon path, is consistently 
found for all optode separations. Thus, the ballistic regime of photon propagation extends to 
about 2 1 . Note, that the uncertainty of the parameter r d increases with the optode distance. 
This shows that the shift ceases to be important as the spacing enters the range of validity of 
the diffusion approximation. The parameter t converges within the given uncertainty to the 
input value l input =4.5 mm. However, the parameter £ a converges to a value which is 
significantly smaller than the input value l a input = 190 mm . This indicates the increasing 
influence of the boundaries on this parameter. Apparently, the optode depth of 40 mm below 
the air interface is too small, causing a cut-off of long photon paths. Thus, the t a measured 
with the real phantom is slightly overestimated (which we take into account in the final error 
estimate in Section 7, Conclusions). However, the resulting error of about 10 mm remains 
well within the experimental error estimate. 

Table 2. Parameters obtained by fitting the diffusion model to the simulated data 
with l" = 4.5 mm, I = 190 mm 

inniil tt mout 



Optode distance 
r [mm] 


Shift r , [ps] 


Transport mean 
path t [mm] 


Absorption length 
l ti [mm] 


24 


50 + 5 


5.5+0.5 


214+16 


36 


50 + 6 


4.8+0.3 


186 + 8 


48 


50+14 


4.6 + 0.2 


179 + 5 


55 


50 + 20 


4.6 + 0.2 


178 + 3 



7. Conclusions 



The optical parameters of our phantom tissue precursor consisting of 1 g L 1 TiC<2 in 16.7% 
gelatin are summarized as follows: 

Refractive index: n lm = 1.359(±0.002) ; 

Transport mean free path: t = 4.5(±0.5) mm [ ju' e * fi' s = (0.22 ± 0.02) mm" 1 ]; 

Absorption length: l a = 190(+10 - 20) mm [ n a = 0.0053(+0.0006 - 0.0003) mm 1 ]; 

Extinction coefficient: fi e » fj, s = 2.86(±0.04) mm 1 . 

The transport mean free path t of almost 5 mm (reduced scattering coefficient ju[ of only 0.2 
mm -1 ) is substantially larger than initially expected for the given concentration of Ti0 2 
particles. This large t reflects the large asymmetry of underlying scattering law. From the 
measured /u[ and jU s we obtain g = 0.93( + 0.01), which is much larger than one would expect 
for Mie scattering from bare Ti0 2 particles (g = 0.5-0.7). Indeed, our preliminary 
measurements with a simple goniometer indicate that the angular dependence of scattering is 
not compatible with the Ti0 2 particles, but can be reasonably modeled with a generalized 
Henyey-Greenstein structure factor (Eq. (9) with g = 0.93 and shape parameter v = -0.3 . We 
speculate that the scattering is due to sparse microstructures in the order of 10 fim in size. 

These findings must be taken into account in the next steps of the phantom development. 
When designing a tissue phantom, one should not blindly expect the optical properties of the 
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product to be the sum of the optical properties of its components. The optical properties 
depend on the microstructure of the system, as it results from interactions between the 
constituents, and whose formation depends on the physico-chemical parameters such as pH 
and ionic strength. Certainly, the microstructure of the gelatin-TiC^ phantoms would deserve a 
thorough investigation. (This, however, would require a dedicated small angle light scattering 
setup.) 

The parameters f and i a (which are the relevant parameters in many practical situations) 
were determined using a novel version of the fiber-optic TOF technique in which we immerse 
the injection and detection optodes deep inside the scattering medium. This geometry 
minimizes the influence of the interfaces on the TOF distribution and allows the analysis of 
the data in terms of the simple free space diffusion approximation. 

The robustness of the approach was tested using our new simulation tool, which allows not 
only realistic modeling of the experimental geometry but also provides a simple way to vary 
the scattering law and its parameters. In the diffusion regime, the parameters obtained from 
simulated and experimentally determined data agree well within their error estimates. The 
simple diffusion approximation can be reliably applied down to optode separations of less 
than 6 I , if one corrects for a small time shift r d « 2f I c l (in our case about 50 ps) that 
accounts for the transition from ballistic to diffusive regime of light propagation. 

Moreover, our simulations confirm that the above four parameters, namely n, f , l a and 
(i e , represent a necessary and sufficient set of optical parameters of a realistic tissue phantom. 
Usually n, t , l a are the most relevant parameters. The fourth parameter fi e (or rather the 
asymmetry parameter g that corresponds to a certain fi e when t is fixed) plays a role only at 
length and time scales that involve a small number of scattering events. (This may be the case, 
for example, with multiphoton imaging of the skin.) In any case, the exact shape of the 
scattering law seems to play a negligible role in most practical situations (although there may 
be exceptions [27]). This certainly simplifies the design and characterization of the phantom 
tissue. 

These insights combined with the experimental and computational tools described here 
create a powerful platform on which to build realistic tissue phantoms. Such optically, 
acoustically, and geometrically realistic phantoms will serve as essential tools in the 
development of important new biomedical imaging techniques. 
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